## Flynn Stewart study 2 (Prolific, UK sample, Nov. 2016)

#setwd("~/Downloads/") #set your working directory here
library(foreign)
data<-read.csv("FlynnStewart study 2 data.csv")
names(data)
dim(data) #n=618

## Recodes:
library(car)

###########################
### Treatment indicator ###
###########################

data$Q188<-recode(data$Q188,"1=1;else=0",as.factor.result=F)
table(data$Q188)
data$Q257<-recode(data$Q257,"1=1;else=0",as.factor.result=F)
table(data$Q257)
data$Q258<-recode(data$Q258,"1=1;else=0",as.factor.result=F)
table(data$Q258)
data$Q259<-recode(data$Q259,"1=1;else=0",as.factor.result=F)
table(data$Q259)
data$Q260<-recode(data$Q260,"1=1;else=0",as.factor.result=F)
table(data$Q260)
data$Q350<-recode(data$Q350,"1=1;else=0",as.factor.result=F)
table(data$Q350)

data$condition<-ifelse(data$Q188==1, 1,
                       ifelse(data$Q257==1, 2,
                              ifelse(data$Q258==1, 3,
                                     ifelse(data$Q259==1, 4,
                                            ifelse(data$Q260==1, 5,
                                                   ifelse(data$Q350==1, 6, NA))))))

table(data$condition) #Ns: 102, 109, 83, 113, 95, 109

# Approval of tactics DV
data$approve.tactics<-recode(data$Q342,"1=7;6=6;2=5;3=4;4=3;5=2;8=1;else=NA",as.factor.result=F)
table(data$approve.tactics)
tapply(data$approve.tactics, data$condition, mean, na.rm=T)

# General approval DV
data$approve.group<-recode(data$Q343,"1=7;2=6;10=5;3=4;6=3;13=2;14=1;else=NA",as.factor.result=F)
table(data$approve.group)
tapply(data$approve.group, data$condition, mean, na.rm=T)

# OLF capable of governing independent state DV
data$capable<-recode(data$Q344, "2=5;5=4;6=3;7=2;8=1;else=NA",as.factor.result=F)
table(data$capable)
tapply(data$capable, data$condition, mean, na.rm=T)

# Support/oppose independent OLF state DV
data$independent<-recode(data$Q345,"1=7;11=6;12=5;13=4;14=3;15=2;16=1;else=NA",as.factor.result=F)
table(data$independent)
tapply(data$independent, data$condition, mean, na.rm=T)

# OLF a legitimate alternative to Ethiopia gov't DV
data$legitimate<-recode(data$Q346, "1=7;2=6;3=5;4=4;5=3;6=2;7=1;else=NA", as.factor.result=F)
table(data$legitimate)
tapply(data$legitimate, data$condition, mean, na.rm=T)

# OLF feeling thermometer DV
data$Q347[1]<-NA
data$Q347<-as.numeric(data$Q347)
data$ft<-recode(data$Q347,"NA=NA",as.factor.result=F)
tapply(data$ft, data$condition, mean, na.rm=T)
summary(data$ft)
hist(data$ft)

# Competence of Ethiopian gov't DV
data$competent<-recode(data$Q319,"1=5;2=4;3=3;4=2;5=1;else=NA", as.factor.result=F)
tapply(data$competent, data$condition, mean, na.rm=T)

# Local residents better/worse off under OLF DV
data$residents<-recode(data$Q320,"1=3;2=1;3=2;else=NA", as.factor.result=F)
tapply(data$residents, data$condition, mean, na.rm=T)

# UN should send peacekeepers DV
data$peacekeepers<-recode(data$Q118_2, "1=7;2=6;3=5;4=4;5=3;6=2;7=1;else=NA", as.factor.result=F)
tapply(data$peacekeepers, data$condition, mean, na.rm=T)

# UK should recognize OLF as rightful gov't DV
data$recognize<-recode(data$Q118_6, "1=7;2=6;3=5;4=4;5=3;6=2;7=1;else=NA", as.factor.result=F)
tapply(data$recognize, data$condition, mean, na.rm=T)

# UK should provide humanitarian assistance to Ethiopian gov't DV
data$assist_govt<-recode(data$Q118_3, "1=7;2=6;3=5;4=4;5=3;6=2;7=1;else=NA", as.factor.result=F)
tapply(data$assist_govt, data$condition, mean, na.rm=T)

# UK should provide humanitarian assistance to OLF DV
data$assist_olf<-recode(data$Q118_4, "1=7;2=6;3=5;4=4;5=3;6=2;7=1;else=NA", as.factor.result=F)
tapply(data$assist_olf, data$condition, mean, na.rm=T)

# UK should provide military/technical support to Ethiopian gov't DV
data$military<-recode(data$Q118_5, "1=7;2=6;3=5;4=4;5=3;6=2;7=1;else=NA", as.factor.result=F)
tapply(data$military, data$condition, mean, na.rm=T)

#########################################
### Covariates / randomization checks ###
#########################################

## Demographics/political:
data$age.temp<-recode(data$Q181,"'In what year were you born?'=NA",as.factor.result=F)
data$age<-(2016 - data$age.temp)
table(data$age)
data$graduated<-recode(data$Q115,"1=0;2=0;3=0;4=1;else=NA",as.factor.result=F)
table(data$graduated)
data$male<-recode(data$Q15,"1=1;else=0",as.factor.result=F)
table(data$male)
data$female<-recode(data$Q15,"2=1;else=0",as.factor.result=F)
table(data$female)
data$white<-recode(data$Q16,"1=1;else=0",as.factor.result=F)
table(data$white)
#PID:
data$pid.3<-recode(data$Q90,"1=1;3=3;4=4;else=NA",as.factor.result=F)
table(data$pid.3)
data$labour<-recode(data$pid.3,"1=1;3:4=0;else=NA",as.factor.result=F)
table(data$labour)
data$conserv<-recode(data$pid.3,"3=1;1=0;4=0;else=NA",as.factor.result=F)
table(data$conserv)
data$libdem<-recode(data$pid.3,"4=1;1=0;3=0;else=NA",as.factor.result=F)
table(data$libdem)
data$ideo.10<-recode(data$Q110,"1=1;2=2;24=3;25=4;26=5;27=6;28=7;29=8;30=9;31=10;else=NA",as.factor.result=F)
table(data$ideo.10)
data$interest<-recode(data$Q54,"1=5;2=4;3=3;4=2;5=1;else=NA",as.factor.result=F)
table(data$interest) 

## General pol. knowledge (count skips as incorrect):
#define high knowledge as above median (5 correct out of 9):

## Rand checks:
tapply(data$age, data$condition, median, na.rm=T)
chisq.test(data$age,data$condition) #p=.19

tapply(data$education, data$condition, median, na.rm=T)
chisq.test(data$education,data$condition) #p=.40

tapply(data$male, data$condition, mean, na.rm=T)
chisq.test(data$male,data$condition) #p=.07

tapply(data$white, data$condition, mean, na.rm=T)
chisq.test(data$white,data$condition) #p=.60

tapply(data$dem, data$condition, mean, na.rm=T)
chisq.test(data$dem,data$condition) #p=.48

tapply(data$rep, data$condition, mean, na.rm=T)
chisq.test(data$rep,data$condition) #p=.68

tapply(data$ideo.7, data$condition, mean, na.rm=T)
chisq.test(data$ideo.7,data$condition) #p=.44


##Regressions for paper (Table 1): 
#create violence and service dummies:
data$violent<-ifelse(data$condition==2 | data$condition==4 | data$condition==6, 1, 0)
table(data$violent, data$condition)
data$restrictive<-ifelse(data$condition==3 | data$condition==4, 1, 0)
table(data$restrictive, data$condition)
data$inclusive<-ifelse(data$condition==5 | data$condition==6, 1, 0)
table(data$inclusive, data$condition)

reg1<-lm(legitimate~violent*restrictive + violent*inclusive, data=data)
summary(reg1)
reg2<-lm(capable~violent*restrictive + violent*inclusive, data=data)
summary(reg2)
reg3<-lm(independent~violent*restrictive + violent*inclusive, data=data)
summary(reg3)

##testing restrictive vs. inclusive in all three models:
linearHypothesis(reg1, "restrictive = inclusive") #n.s.
linearHypothesis(reg2, "restrictive = inclusive") #n.s.
linearHypothesis(reg3, "restrictive = inclusive") #n.s.

##testing treatment effects on dichotomized outcomes (legitimate/not legitimate:) - for appendix
data$legitimate.b<-ifelse(data$legitimate %in% c(5,6,7),1,0)
table(data$legitimate,data$legitimate.b)
data$capable.b<-ifelse(data$capable %in% c(4,5),1,0)
table(data$capable,data$capable.b)
data$independent.b<-ifelse(data$independent %in% c(5,6,7),1,0)
table(data$independent,data$independent.b)

reg4<-lm(legitimate.b~violent*restrictive + violent*inclusive, data=data)
summary(reg4)
reg5<-lm(capable.b~violent*restrictive + violent*inclusive, data=data)
summary(reg5)
reg6<-lm(independent.b~violent*restrictive + violent*inclusive, data=data)
summary(reg6)

reg7<-glm(legitimate.b~violent*restrictive + violent*inclusive, data=data, family=binomial(link="logit"))
summary(reg7)
reg8<-glm(capable.b~violent*restrictive + violent*inclusive, data=data, family=binomial(link="logit"))
summary(reg8)
reg9<-glm(independent.b~violent*restrictive + violent*inclusive, data=data, family=binomial(link="logit"))
summary(reg9)


##Confidence intervals for Violence x Services interaction effects:
#NOTE: confidence intervals around these effects were used to make error bars in Figure 1
#make two treatments factors:
data$violence.fac<-ifelse(data$condition==1 | data$condition==3 | data$condition==5, "base-nv", "v")
data$violence.fac<-as.factor(data$violence.fac)
data$service.fac<-ifelse(data$condition==1 | data$condition==2, "base-none",
                         ifelse(data$condition==3 | data$condition==4, "rest", 
                                ifelse(data$condition==5 | data$condition==6, "incl", NA)))
data$service.fac<-as.factor(data$service.fac)
table(data$condition, data$violence.fac) #confirm 
table(data$condition, data$service.fac) #confirm
#Marginal effects of (violence | services) on LEGITIMACY DV:
reg1<-lm(legitimate~violence.fac*service.fac, data=data)
summary(reg1) #coefs for Table in paper
library(contrast) #calculate marginal effects using "contrast" package:
#Effect of violence among groups that provide NO services:
ViolenceEffect_NoServices <- contrast(reg1,
                                      list(service.fac = "base-none", violence.fac = "v"),
                                      list(service.fac = "base-none", violence.fac = "base-nv"))
print(ViolenceEffect_NoServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide RESTRICTIVE services:
ViolenceEffect_RestServices <- contrast(reg1,
                                        list(service.fac = "rest", violence.fac = "v"),
                                        list(service.fac = "rest", violence.fac = "base-nv"))
print(ViolenceEffect_RestServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide INCLUSIVE services:
ViolenceEffect_InclServices <- contrast(reg1,
                                        list(service.fac = "incl", violence.fac = "v"),
                                        list(service.fac = "incl", violence.fac = "base-nv"))
print(ViolenceEffect_InclServices, X = T) #Contrast=effect of violence

#looking at extent of attenuation in violence effect across service provision 
#(t-test: violent no services vs. violent inclusive services):
t.test(data$legitimate[data$violent==1 & data$restrictive==0 & data$inclusive==0],
       data$legitimate[data$violent==1 & data$restrictive==0 & data$inclusive==1], na.rm=T)

#Marginal effects of (violence | services) on CAPABLE OF GOVERNING DV:
reg2<-lm(capable~violence.fac*service.fac, data=data)
summary(reg2) #coefs for Table 4 in paper
#Effect of violence among groups that provide NO services:
ViolenceEffect_NoServices <- contrast(reg2,
                                      list(service.fac = "base-none", violence.fac = "v"),
                                      list(service.fac = "base-none", violence.fac = "base-nv"))
print(ViolenceEffect_NoServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide RESTRICTIVE services:
ViolenceEffect_RestServices <- contrast(reg2,
                                        list(service.fac = "rest", violence.fac = "v"),
                                        list(service.fac = "rest", violence.fac = "base-nv"))
print(ViolenceEffect_RestServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide INCLUSIVE services:
ViolenceEffect_InclServices <- contrast(reg2,
                                        list(service.fac = "incl", violence.fac = "v"),
                                        list(service.fac = "incl", violence.fac = "base-nv"))
print(ViolenceEffect_InclServices, X = T) #Contrast=effect of violence

#Marginal effects of (violence | services) on SUPPORT INDEPENDENT STATE DV:
reg3<-lm(independent~violence.fac*service.fac, data=data)
summary(reg3) #coefs for Table 4 in paper
#Effect of violence among groups that provide NO services:
ViolenceEffect_NoServices <- contrast(reg3,
                                      list(service.fac = "base-none", violence.fac = "v"),
                                      list(service.fac = "base-none", violence.fac = "base-nv"))
print(ViolenceEffect_NoServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide RESTRICTIVE services:
ViolenceEffect_RestServices <- contrast(reg3,
                                        list(service.fac = "rest", violence.fac = "v"),
                                        list(service.fac = "rest", violence.fac = "base-nv"))
print(ViolenceEffect_RestServices, X = T) #Contrast=effect of violence 
#Effect of violence among groups that provide INCLUSIVE services:
ViolenceEffect_InclServices <- contrast(reg3,
                                        list(service.fac = "incl", violence.fac = "v"),
                                        list(service.fac = "incl", violence.fac = "base-nv"))
print(ViolenceEffect_InclServices, X = T) #Contrast=effect of violence


##Export regression tables:
library(stargazer)
stargazer(reg1, reg2, reg3, type="latex",  title="OLS Regression Models Predicting Legitimacy Evaluations (Study 1)",
          covariate.labels=c("Violent", "Restrictive Services","Inclusive Services", "Violent x Restrictive Services","Violent x Inclusive Services","Constant"), omit.stat=c("rsq", "adj.rsq", "res.dev", "F"),
          dep.var.labels = c("Legitimate Alternative","Capable of Governing","Support Independent State"),
          model.numbers=F, digits=2)

#robustness check: dichotomized outcomes with OLS (for appendix):
stargazer(reg1, reg2, reg3, reg4, reg5, reg6, type="latex",  title="Regression Models Predicting Legitimacy Evaluations (Study 1)",
          covariate.labels=c("Violent", "Restrictive Services","Inclusive Services", "Violent x Restrictive Services","Violent x Inclusive Services","Constant"), omit.stat=c("rsq", "adj.rsq", "res.dev", "F"),
          dep.var.labels = c("Legitimate Alternative (1-7)","Capable of Governing (1-5)","Support Independent State (1-7)",
                             "Legitimate Alternative (0-1)", "Capable of Governing (0-1)", "Support Independent State (0-1)"),
          model.numbers=F, digits=2)
